#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <algorithm>
#include <ctime>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <mpi.h>
#include <unistd.h>
#include <math.h>
#include <exception>

#include "PNN1.h"

#define SM (CLS / sizeof (double))

int confusion[30][30];
int confusion_sum[30][30];
int contactStat[30];
int contactStat_sum[30];
int sum[30];
int test_mode=0;

SEQUENCE::SEQUENCE(int length_seq, PNN1* pModel, string name)
{
  m_pModel = pModel;
  p_name= name;
  this->length_seq = length_seq;
  if (test_mode==0)
    length_contact=length_seq*length_seq/1;
  else
    length_contact=length_seq*length_seq;
  dim_features=m_pModel->dim_features;
  window_size= m_pModel->window_size;
  dim_one_pos=m_pModel->dim_one_pos;
  num_gates = m_pModel->num_gates;
  num_states = m_pModel->num_states;
  num_labels = m_pModel->num_labels;
  dim_pairwise_features = m_pModel->dim_pairwise_features;
  gates=NULL;
  arrVi=NULL;
  Partition=NULL;
  obs_feature = new Score*[length_seq];
  for(int i=0;i<length_seq;i++)
    obs_feature[i] = new Score[dim_one_pos];
  s_start=0;
  s_end=length_seq-1;
  output_center=NULL;
  predicted_label=NULL;
  obj_weights=NULL;
  Rg=2.2*pow(length_seq, 0.38);
  pFeatures=new Score[dim_features];
}

SEQUENCE::~SEQUENCE(){
  delete arrVi;
  delete gates;
  delete Partition;
  delete obs_label;

  for(int i=0;i<length_seq;i++)
    delete obs_feature[i];
  delete obs_feature;
  delete output_center;
  delete predicted_label;
  delete pFeatures;
  delete obj_weights;
}

Score SEQUENCE::Obj(){ // observed label : log-likelihood
  ComputeGates();
  ComputeVi();
  CalcPartition();

  Score obj = 0;
  for(int t=0;t<nDataPoints;t++)
    obj+= obj_weights[t] * ( (*arrVi)(t, obs_label[t]) - Partition[t] );
  
  delete arrVi; arrVi=NULL;
  //if (proc_id==0) cerr << "obj=" << obj << "; Partition=" << Partition << endl;
  return obj;
}

Score SEQUENCE::MAP(){ // Posterior Decoding (Marginal Probability Decoder)
  //memset(pred, 0, sizeof(int)*20);
  Score bestScore, obj = 0;
  Score *vi2;
  int i, t, bestState=0; 

  CalcPartition();
  for(t=nDataPoints-1; t>=0;t--)
    obj+= obj_weights[t] * ( (*arrVi)(t, obs_label[t]) - Partition[t] );

  for(t=nDataPoints-1, vi2=&((*arrVi)(t, 0)); t>=0; t--, vi2-=arrVi->cols){
    bestState=0; 
    bestScore=vi2[0];
    for(i=1; i< num_states; i++){
      if (vi2[i]>bestScore){
        bestScore=vi2[i];
        bestState=i;
      }
    }
    predicted_label[t]=bestState;
    //pred[predicted_label[t]]++;
  }
  return obj;
}

void SEQUENCE::CalcPartition()
{
  int t, i, t2, i2, lastt=nDataPoints/SM*SM, lasti=num_states/SM*SM;
  Score score, *vi2, *P2;

  if (Partition==NULL) 
    Partition = new Score[nDataPoints];
  for(t=0;t<nDataPoints;t++)
    Partition[t] = (Score)LogScore_ZERO;
// LogScore_PLUS_EQUALS(Partition[t], (*arrVi)(t, state));
  for(t=0;t<nDataPoints;t+=SM)
    for(i=0; i< num_states; i+=SM)
      for ( t2 = 0, vi2 = &((*arrVi)(t, i)), P2 = &Partition[t];
            t2 < (t<lastt?SM:nDataPoints%SM); t2++, vi2 += arrVi->cols)
        for (i2=0; i2 < (i<lasti ? SM : num_states%SM); i2++){
          score = vi2[i2];
          LogScore_PLUS_EQUALS(P2[t2], score);
          //if (proc_id==0) cerr << "  Partition(" << t << "," << state << ")=" << score << " : " << Partition[t] << endl;
        }
  //cerr << endl;
}

// Calculate phi_\theta(x_i, y)
void SEQUENCE::ComputeVi()
{
  int pos, pos2, l, l2, g, g2, 
      lastpos=nDataPoints/SM*SM, lastl=num_states/SM*SM, lastg=num_gates/SM*SM;
  Score *gates2, *vi2, *pW2;

  if (arrVi==NULL)
    arrVi = new ScoreMatrix(nDataPoints, num_states, "arrVi");
  arrVi->resize(nDataPoints, num_states); //arrVi->Fill(LogScore_ZERO);

// (*arrVi)(pos,l) += m_pModel->weights[GetWeightLAddr(l, g)]*(*gates)(pos,g);

  for ( pos=0; pos<nDataPoints; pos+=SM)
  for ( l=0; l<num_states; l+=SM)
  for ( g=0; g<num_gates; g+=SM) {
    for ( pos2 = 0, gates2=&((*gates)(pos, g)), vi2=&((*arrVi)(pos, l)); 
          pos2 < (pos<lastpos?SM:nDataPoints%SM);
          pos2++, vi2+=arrVi->cols, gates2+=gates->cols)
      for ( l2 = 0, pW2=m_pModel->GetWeightLAddr(l,g); 
            l2 < (l<lastl ? SM : num_states%SM); 
            l2++, pW2 += num_gates){
        for ( g2 = 0; g2 < (g<lastg ? SM : num_gates%SM); g2++) {
          vi2[l2] += pW2[g2] * gates2[g2];
          // if (vi2[l2]>10000)cerr<<" ComputeVi: Vi("<<pos+pos2<<", "<<l+l2<<")= W("<<l+l2<<","<<g+g2<<") * gates("<<pos+pos2<<","<<g+g2<<"):: "<<vi2[l2]<<" += "<<pW2[g2]<<" * "<<gates2[g2]<<endl; 
        } // exit(0);
      }
  }
  
  //for (pos=0; pos<length_contact; pos++)
  //  for (l=0; l<num_states; l++)
  //    cerr << "  ComputeVi(" << pos << ", " << l << ")=" << (*arrVi)( pos, l) << endl;
}

SingletonSequence::SingletonSequence(int len, PNN1* pModel, string name)
  : SEQUENCE(len, pModel, name)
{
  obs_label = new int[length_seq];
  _features=NULL;
}

SingletonSequence::~SingletonSequence()
{
  delete _features;
}

int SingletonSequence::makeFeatures(){
  //build features for local windows  
  _features = new ScoreMatrix(length_seq, dim_features, "_features");
  for(int t=0;t<length_seq;t++){
    if (t<window_size/2)
      (*_features)(t, t) = 1;
    else if (length_seq-1-t<window_size/2)
      (*_features)(t, length_seq-1-t) = 1;

    int pivot = window_size/2*2;
    for(int i=0; i< window_size; i++){
      int offset = t + i - window_size/2;
      if(offset <0 || offset >=length_seq){
        for(int j=0; j<dim_one_pos; j++)
          (*_features)(t, pivot) = 0, pivot++;
      } else {
        for(int j=0; j<dim_one_pos; j++)
          (*_features)(t, pivot) = obs_feature[offset][j], pivot++;
      }
    }
    (*_features)(t, pivot) = 1; // last one for the bias term
    // cerr << " ### last pivot=" << pivot << endl;
  }

  nDataPoints = length_seq;
  predicted_label = new int[nDataPoints];
  obj_weights=new double[nDataPoints];
  for (int t=0; t<nDataPoints; t++)
    obj_weights[t]=1.0;

  if (proc_id==0) {
    for (int i =0; i<dim_features; i++) cerr << i << ":" << (*_features)(0, i) << " ";
    cerr << endl;
  }
  return 0;
}

void SingletonSequence::ComputeGates()
{
  Score *pfb2, *ocenter, *pW2;
  int t, t2, g, g2, k, k2,
      lasttseq=length_seq/SM*SM, lastg=num_gates/SM*SM, lastk=dim_features/SM*SM;

  struct timeval ts, te;
  gettimeofday(&ts, 0);

  if (gates==NULL){
    gates = new ScoreMatrix(length_seq, num_gates, "gates");
    output_center = new ScoreMatrix(num_gates,length_seq, "output_center");
  } else {
    output_center->resize(num_gates, length_seq); //output_center->Fill(0);
  }
  for (t=0; t<length_seq; t+=SM)
  for (g=0; g<num_gates; g+=SM)
  for (k=0; k<dim_features; k+=SM) 
    for ( g2 = 0, pW2=m_pModel->GetWeightGAddr(g,k), ocenter=&((*output_center)(g, t)); 
          g2 < (g<lastg?SM:num_gates%SM); 
          g2++, pW2+=dim_features, ocenter+=output_center->cols) {
      for ( t2 = 0, pfb2=&((*_features)(t, k)); 
            t2<(t<lasttseq?SM:length_seq%SM); 
            t2++, pfb2+=_features->cols) 
      {
        for (k2=0; k2<(k<lastk?SM:dim_features%SM); k2++) {
          ocenter[t2] += pW2[k2] * pfb2[k2];
        }
      } //cerr << "  GetGate(" << gate << ", " << pos << ")=" << output << endl; //exit(0);
    }
  for (t=0; t<length_seq; t++)
    for (g=0; g<num_gates; g++){
      (*gates)(t, g)=Gate((*output_center)(g, t)); // if (t==30) cerr << "   gate("<< t <<"," << g <<")= " << (*gates)(t, g) << endl;
    }
if (proc_id==0) {
  gettimeofday(&te, 0);
  double dif = difftime (te.tv_sec,ts.tv_sec);
  //cerr << "Gates: " << dif*1000+te.tv_usec/1000-ts.tv_usec/1000 << " miliseconds." << endl;
} //exit(0);
}

Score SingletonSequence::Predict()
{
  ComputeGates();
  ComputeVi();
  Score obj =MAP();
  ComputeTestAccuracy();
  return obj;
}

Score SingletonSequence::ComputeGradient(bool bCalculateGate)
{
  int gstart = m_pModel->iGateStart;
  ScoreMatrix* prob_weight_sum=new ScoreMatrix(length_seq, num_gates, "prob_weight_sum");
  ScoreMatrix* probs=new ScoreMatrix(length_seq, num_states, "probs");
  double *pGrad2, *pGate2, *pW2, *probs2;
  Score *pPWS2, out, prob_out;
  int weightIdx, t, t2, g, g2, l, l2, k, k2, 
      lastt=length_seq/SM*SM, lastg=num_gates/SM*SM, 
      lastl=num_states/SM*SM, lastk=dim_features/SM*SM;

  ScoreMatrix* delta_buffer = new ScoreMatrix(num_gates, length_seq, "delta_buffer");
  Score *pdb2, *pfb2;

  // Compute Gate, Vi and partition function;
  struct timeval ts, te;
  gettimeofday(&ts, 0);

  ComputeGates();
  ComputeVi();
  CalcPartition();
/////////////////////////////////Test for Online Learning /////////////////
  Score s, obj = 0;
  for(t=0;t<nDataPoints;t++){
    s=(*arrVi)(t, obs_label[t]);
    obj+=s-Partition[t];
  }
////////////////////////////////////////////////////////////////////////////

if (proc_id==0) {
  gettimeofday(&te, 0);
  double dif = difftime (te.tv_sec,ts.tv_sec);
  //cerr << "ComputeGateViPartition: " << dif*1000+te.tv_usec/1000-ts.tv_usec/1000 << " miliseconds." << endl;
  gettimeofday(&ts, 0);
}
  // compute the empirical values
  prob_weight_sum->Fill(0);
  for(t=0; t<length_seq; t++){
    weightIdx=gstart + obs_label[t]*num_gates;
    for(g=0, pW2=&m_pModel->weights[weightIdx], pGrad2=&m_pModel->grad[weightIdx]; g<num_gates; g++){
      out = (*gates)(t, g);
      pGrad2[g] += out; // label-based gates weights
      (*prob_weight_sum)(t,g) += pW2[g]*(1.0-out)*out; // Gate weights; h'(x)=[1-h(x)]*h(x)=(1.0-out)*out
    }
  }
  for(g=0;g<num_gates;g++)
    for(t=0;t<nDataPoints;t++)
      (*delta_buffer)(g, t)+= (*prob_weight_sum)(t,g);
  for(g=0; g<num_gates; g+=SM)
  for(k=0; k<dim_features; k+=SM)
  for(t=0; t<length_seq; t+=SM)
    for ( g2 = 0, pdb2=&((*delta_buffer)(g,t)), 
                  pGrad2 = &m_pModel->grad[m_pModel->iFeatureStart + g*dim_features+k]; 
          g2 < (g<lastg?SM:num_gates%SM); 
          g2++, pdb2+=delta_buffer->cols, pGrad2+=dim_features)
    for ( t2 = 0, pfb2 = &((*_features)(t, k)); 
          t2 < (t<lastt ? SM : length_seq%SM); 
          t2++, pfb2+=_features->cols) 
    {
      for ( k2=0; k2<(k<lastk?SM:dim_features%SM); k2++)
      {
        pGrad2[k2] += pdb2[t2]*pfb2[k2];
      }
    }

if (proc_id==0) {
  gettimeofday(&te, 0);
  double dif = difftime (te.tv_sec,ts.tv_sec);
  //cerr << "Gradient Empirical: " << dif*1000+te.tv_usec/1000-ts.tv_usec/1000 << " miliseconds." << endl;
  gettimeofday(&ts, 0);
}
  // compute the expected values
  prob_weight_sum->Fill(0);
  for(t=0;t<length_seq;t++)
    for (l=0; l<num_states; l++)
      (*probs)(t,l) = exp( (*arrVi)(t, l) - Partition[t] );

// out=(*gates)(t, g); prob_out = prob(t, l)*out; 
// m_pModel->grad[l, g] -= prob_out; 
// (*prob_weight_sum)(t,g) += prob_out*weightsL(l, g)*(1.0-out);

  for ( t=0;t<length_seq;t+=SM)
  for ( l=0; l<num_states; l+=SM)
  for ( g=0; g<num_gates; g+=SM)
    for ( t2 = 0, probs2=&((*probs)(t,l)), 
                  pPWS2=&((*prob_weight_sum)(t,g)), 
                  pGate2=&((*gates)(t,g)); 
          t2 < (t<lastt ? SM : length_seq%SM); 
          t2++, probs2+=probs->cols, 
                pPWS2+=prob_weight_sum->cols, 
                pGate2+=gates->cols) 
      for ( l2 = 0, pW2 = m_pModel->GetWeightLAddr(l, g),
                    pGrad2 = &m_pModel->grad[gstart+l*num_gates+g];
            l2 < ( l<lastl ? SM : num_states%SM); 
            l2++, pW2 += num_gates, pGrad2 += num_gates){
        for ( g2 = 0; g2 < (g<lastg ? SM : num_gates%SM); g2++){
          prob_out = probs2[l2] * pGate2[g2];
          pGrad2[g2] -= prob_out;
          pPWS2[g2]+= prob_out * pW2[g2] * (1.0-pGate2[g2]);
        }
      }

  delta_buffer->Fill(0);
  for(g=0;g<num_gates;g++){
    for(t=0;t<length_seq;t++){
      (*delta_buffer)(g, t)+= (*prob_weight_sum)(t,g);
    }
  }
  for(g=0; g<num_gates; g+=SM)
  for(k=0; k<dim_features; k+=SM)
  for(t=0; t<length_seq; t+=SM)
    for ( g2 = 0, pdb2 = &((*delta_buffer)(g,t)), 
                  pGrad2 = &m_pModel->grad[m_pModel->iFeatureStart + g*dim_features + k];
          g2 < (g<lastg ? SM : num_gates % SM);
          g2++, pdb2 += delta_buffer->cols, pGrad2 += dim_features)
    for ( t2 = 0, pfb2 = &((*_features)(t, k));
          t2 < (t<lastt ? SM : length_seq%SM); 
          t2++, pfb2 += _features->cols) 
    {
      for (k2 = 0; k2 < (k<lastk ? SM : dim_features%SM); k2++){
        pGrad2[k2] -= pdb2[t2] * pfb2[k2];
      }
    }

if (proc_id==0) {
  gettimeofday(&te, 0);
  double dif = difftime (te.tv_sec,ts.tv_sec);
  //cerr << "Gradient Expected: " << dif*1000+te.tv_usec/1000-ts.tv_usec/1000 << " miliseconds." << endl;
}

  delete delta_buffer;
  delete probs,
  delete prob_weight_sum;
  delete arrVi; arrVi=NULL;

  return obj;
}

void SingletonSequence::ComputeTestAccuracy(){
  m_pModel->totalPos += nDataPoints;
  // comparison
  for(int t=0; t<nDataPoints;t++) {
    if(obs_label[t]==predicted_label[t])
      m_pModel->totalCorrect++;
    confusion[predicted_label[t]][obs_label[t]]++;
  }
  //for(int t=0; t<nDataPoints;t++) cerr << predicted_label[t] << ":" << obs_label[t] << " "; 
  //cerr << "=="<<m_pModel->totalCorrect << endl;
}

PairwiseSequence::PairwiseSequence(int len, PNN1* pModel, string name)
  : SEQUENCE(len, pModel, name)
{
  contact_DIST=7;
  IdxMap=new int[length_seq*length_seq];
  memset(IdxMap, -1, sizeof(int)*length_seq*length_seq);
  rIdxMap=new int[length_seq*length_seq];
  centerMap=new int[length_seq*length_seq];
  IdxMapCenter=new int[length_seq];
  IdxMapCofactor=new int[length_seq];
  for (int i=0; i<length_seq; i++)
    IdxMapCofactor[i]=IdxMapCenter[i]=-1;
  output_cofactor=NULL;
  distances=new ScoreMatrix(length_seq, length_seq, "distances");
  distances->Fill(-1);
  obs_label=NULL;
  pairwise_features=NULL;
  obs_pairwise_features=NULL;
  if (dim_pairwise_features)
    obs_pairwise_features=new Array3D<double>(length_seq, length_seq, dim_pairwise_features);
}

PairwiseSequence::~PairwiseSequence()
{
  delete output_cofactor;

  delete IdxMap; 
  delete rIdxMap;
  delete centerMap;
  delete IdxMapCenter;
  delete IdxMapCofactor;
  delete distances;
  delete obs_pairwise_features;
  delete pairwise_features;
}

int PairwiseSequence::makeFeatures(){
  if (obs_label==NULL) 
    obs_label = new int[length_contact];
  int df = length_contact*2;
  //cerr << proc_id << "@ df=" << df << ", len_contact=" << length_contact << ", gaps_size=" << gaps.size();
  _features = new int[df];
  
  int totalContact=0;

  //cerr << " window_size=" << window_size << ", len_seq=" << length_seq << endl;
  //build features for local windows  
  for(int t=0;t<length_seq-contact_DIST;t++){
      bool bGap=false;
      for (int g=0; g<(int)gaps.size(); g++){
        if (abs(gaps[g]-t)<=window_size/2){
      bGap=true; break;
    }
      }
      if(bGap)
    continue;
      //for(int u=t+7;u<length_seq-window_size/2;u++){
      for(int u=t+contact_DIST;u<length_seq;u++){
    if (totalContact>=length_contact-1)
      break;
    bGap=false;
    for (int g=0; g<(int)gaps.size(); g++){
      if (abs(gaps[g]-u)<=window_size/2){
        bGap=true; break;
      }
    }
    if(bGap) 
      continue;
    int label=(int)(*distances)(t,u);
    if (label<0)
      continue;
    obs_label[totalContact]= label;
    //if (!proc_id)
    //  cerr << "    Adding contact " << totalContact << ":" << p_name 
    //       << " (" << t << ", " << u << ") with dist=" << (*distances)(t,u)
    //       << " labeled " << obs_label[totalContact] << endl;
    centerMap[totalContact]=t;
    IdxMapCenter[t]=totalContact;
    rIdxMap[totalContact]=u;
    IdxMapCofactor[u]=totalContact;
    int pivot = totalContact*2;
    totalContact++;
    int offset1 = t-window_size/2;
    int offset2 = u-window_size/2;
    _features[pivot++]=offset1;
    _features[pivot++]=offset2;
      }
  }
  delete distances; distances=NULL;
  nDataPoints=length_contact=totalContact; 
  if (length_contact==0)
     return -1; //cerr << proc_id << " Making features..." << p_name<< endl;
  predicted_label = new int[nDataPoints];
  pairwise_features = new ScoreMatrix(dim_pairwise_features+1, nDataPoints, "pairwise_features");
  for(int t=0;t<length_contact; t++){
    int offset1=centerMap[t], offset2=rIdxMap[t];
    (*pairwise_features)(0,t)= (offset2-offset1+0.000001)/length_seq;
    for (int d=1; d<=dim_pairwise_features; d++)
      (*pairwise_features)(d,t)=(*obs_pairwise_features)(offset1, offset2, d-1);
  }
  delete obs_pairwise_features;
  obs_pairwise_features = NULL;
  obj_weights=new double[nDataPoints];
  for (int t=0; t<nDataPoints; t++)
    obj_weights[t]=1.0;
  return 0;
}

// Features: ends vector(2*[window_size/2]), features read from input 
// file (dim_one_pos*window_size), sequence distance; sequence distance/length_seq.
Score PairwiseSequence::getFeatures(int pos, int i){
  int pivot=pos*2;
  int offset1 = _features[pivot];
  int offset2 = _features[pivot+1];
  int iif=i-window_size/2*2;
  int idx = iif / m_pModel->dim_one_pos;
  Score ret=1; // when i==dim_features-3-dim_pairwise_features
  if (i<window_size/2*2){
    offset1+=window_size/2;
    offset2+=window_size/2;
    //if (proc_id==0) cerr << "offset1="<<offset1<<", offset2="<<offset2<<endl;
    if (offset1<window_size/2 && i==offset1)
      return 1;
    else if (length_seq-offset2==window_size/2*2-i && offset2>=length_seq-window_size/2){
      //if (proc_id==0) cerr << length_seq << "-" << offset2 << "="<<window_size/2*2<<"-" <<i << endl; 
      return 1;
    } else 
      return 0;
  } else if (idx<window_size){
    if (offset1+idx<0 || offset1+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset1+idx][iif%m_pModel->dim_one_pos];
  } else if (i<dim_features-3-dim_pairwise_features){
    idx -= m_pModel->window_size;
    // cerr << idx << " - " << i << "*"<<dim_features <<" : " 
    //      << i%m_pModel->dim_one_pos<<" : " << length_seq 
    //      << " : " << offset1 << "  : " << offset2;
    if (offset2+idx<0 || offset2+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset2+idx][iif%m_pModel->dim_one_pos];
    //cerr << " : " << ret << endl;
  } else if (i==dim_features-2-dim_pairwise_features)
    ret = Rg/10;
  else if (i==dim_features-1-dim_pairwise_features)
    ret = (offset2-offset1+0.000001)/length_seq;

  return ret;
}

void PairwiseSequence::ComputeGates()
{
  ScoreMatrix* feature_buffer = new ScoreMatrix(dim_features, length_seq, "feature_buffer");
  ScoreMatrix* feature_buffer_cofactor = new ScoreMatrix(dim_features, length_seq, "feature_buffer_cofactor");
  Score *pfb2, *pfbc2, *ocenter, *ocofactor;
  //int pos;
  int gate, p, t, t2, g, g2, k, k2, kk, 
      lastt = length_contact/SM*SM, lasttseq = length_seq/SM*SM, 
      lastg = num_gates/SM*SM, lastk = (dim_features-1-dim_pairwise_features)/SM*SM;
  Score* gates2, *pW2, output;

  struct timeval ts, te;
  gettimeofday(&ts, 0);

  if (gates==NULL){
    gates = new ScoreMatrix(length_contact, num_gates, "gates");
    output_center = new ScoreMatrix(num_gates,length_seq, "output_center");
    output_cofactor = new ScoreMatrix(num_gates,length_seq, "output_cofactor");
  } else {
    // gates->resize(length_contact, num_gates); //gates->Fill(LogScore_ZERO);
    output_center->resize(num_gates, length_seq); //output_center->Fill(0);
    output_cofactor->resize(num_gates, length_seq); // output_cofactor->Fill(0);
  }
  for(t=0; t<length_seq; t++){
    int idx=IdxMapCenter[t];
    int coidx=IdxMapCofactor[t];
    for(k=0;k<dim_features-1;k++){
      if (k<window_size/2 || (k>=window_size/2*2 && k<dim_features/2+window_size/2-2)){
        if (idx<0) { continue;}
        (*feature_buffer)(k, t)= getFeatures(idx, k);
      } else {
        if (coidx<0) continue;
        (*feature_buffer_cofactor)(k, t) = getFeatures(coidx, k);
      }
    }
  }
  for (t=0; t<length_seq; t+=SM) 
  for (g=0; g<num_gates; g+=SM)
  for (k=0; k<dim_features-1; k+=SM) 
    for (g2=0, pW2=m_pModel->GetWeightGAddr(g,k), 
               ocenter=&((*output_center)(g, t)), 
               ocofactor=&((*output_cofactor)(g, t)); 
         g2<(g<lastg?SM:num_gates%SM); 
         g2++) 
    {
      for ( k2 = 0, pfb2 = &((*feature_buffer)(k, t)),
                    pfbc2 = &((*feature_buffer_cofactor)(k, t)); 
            k2 < (k<lastk ? SM : (dim_features-1-dim_pairwise_features)%SM); 
            k2++) 
      {
        kk = k + k2;
        if (kk<window_size/2 || (kk>=window_size/2*2 && kk<dim_features/2+window_size/2-2)){
          for(t2=0; t2<(t<lasttseq?SM:length_seq%SM); t2++) 
            ocenter[t2] += pW2[k2] * pfb2[t2];
        } else {
          for(t2=0; t2<(t<lasttseq?SM:length_seq%SM); t2++)
            ocofactor[t2] +=pW2[k2] * pfbc2[t2];
        }
        pfb2+=feature_buffer->cols;
        pfbc2+=feature_buffer_cofactor->cols;
      } //cerr << "  GetGate(" << gate << ", " << pos << ")=" << output << endl; //exit(0);
      pW2+=dim_features;
      ocenter+=output_center->cols;
      ocofactor+=output_cofactor->cols;
    }
  for (t=0; t<length_contact; t+=SM)
  for (g=0; g<num_gates; g+=SM)
    for (t2=0, gates2=&((*gates)(t, g)); 
         t2<(t<lastt?SM:length_contact%SM); 
         t2++, gates2+=gates->cols)
    {
      int center=centerMap[t+t2];
      int cofactor=rIdxMap[t+t2];
      for (g2=0; g2<(g<lastg?SM:num_gates%SM); g2++) {
        gate=g+g2; //pos=t+t2;
        pW2=m_pModel->GetWeightGAddr(gate,0);

        output=(*output_center)(gate, center) + (*output_cofactor)(gate, cofactor);
        for (p=0; p<dim_pairwise_features+1; p++)
          output += pW2[dim_features-dim_pairwise_features-1+p] * (*pairwise_features)(p,t);

        gates2[g2] = Gate(output);
        //if (pos=30) cerr << "   gate("<< gate <<"," << pos <<")= " << gates2[g2] << " " << output << endl;
      }
    }
if (proc_id==0) {
  gettimeofday(&te, 0);
  double dif = difftime (te.tv_sec,ts.tv_sec);
  //cerr << "Gates: " << dif*1000+te.tv_usec/1000-ts.tv_usec/1000 << " miliseconds." << endl;
} //exit(0);
  delete feature_buffer;
  delete feature_buffer_cofactor;
}

Score PairwiseSequence::Predict()
{
  ComputeGates();
  ComputeVi();
  Score obj =MAP();
  ComputeTestAccuracy();
  return obj;
}

Score PairwiseSequence::ComputeGradient(bool bCalculateGate)
{
  int gstart = m_pModel->iGateStart;
  int window_size=m_pModel->window_size;
  ScoreMatrix* prob_weight_sum = new ScoreMatrix(length_contact, 
                                                 num_gates, "prob_weight_sum");
  ScoreMatrix* probs=new ScoreMatrix(length_contact, num_states, "probs");
  double *pGrad2, *pGate2, *pW2, *probs2;
  Score *pPWS2, out, prob_out;
  int weightIdx, p, t, t2, g, g2, l, l2, k, k2, kk,
      lasttseq = length_seq/SM*SM, lastt = length_contact / SM * SM, 
      lastg = num_gates/SM*SM, lastl = num_states / SM * SM,
      lastk = (dim_features-1-dim_pairwise_features) / SM * SM;

  ScoreMatrix* delta_buffer = new ScoreMatrix(num_gates, 
                                              length_seq, "delta_buffer");
  ScoreMatrix* delta_buffer_cofactor = new ScoreMatrix(num_gates,
                                              length_seq, "delta_buffer_cofactor");
  ScoreMatrix* feature_buffer = new ScoreMatrix(dim_features, 
                                              length_seq, "feature_buffer");
  ScoreMatrix* feature_buffer_cofactor = new ScoreMatrix(dim_features,
                                              length_seq, "feature_buffer_cofactor");
  Score *pdb2, *pdbc2, *pfb2, *pfbc2;

  // Compute Gate, Vi and partition function;
  struct timeval ts, te;
  gettimeofday(&ts, 0);
  //if (bCalculateGate) {
    //cerr << " grad--Gates" << endl;
    ComputeGates();
    ComputeVi();
  //}

  CalcPartition();
/////////////////////////////////Test for Online Learning /////////////////
  Score s, obj = 0;
  for(t=0;t<length_contact;t++){
    s=(*arrVi)(t, obs_label[t]);
    obj+=s-Partition[t];
  }
////////////////////////////////////////////////////////////////////////////

if (proc_id==0) {
  gettimeofday(&te, 0);
  double dif = difftime (te.tv_sec,ts.tv_sec);
  //cerr << "ComputeGateViPartition: " << dif*1000+te.tv_usec/1000-ts.tv_usec/1000 << " miliseconds." << endl;
  gettimeofday(&ts, 0);
}
  // compute the empirical values
  prob_weight_sum->Fill(0);
  for(t=0; t<length_contact; t++){
    weightIdx=gstart + obs_label[t]*num_gates;
    for ( g = 0, pW2 = &m_pModel->weights[weightIdx], pGrad2 = &m_pModel->grad[weightIdx]; 
          g < num_gates; g++)
    {
      out = (*gates)(t, g);
      pGrad2[g] += obj_weights[t]*out; // label-based gates weights
      (*prob_weight_sum)(t,g) += obj_weights[t]*pW2[g]*(1.0-out)*out; // Gate weights; h'(x)=[1-h(x)]*h(x)=(1.0-out)*out
    }
  }
  for(g=0;g<num_gates;g++){
    for(t=0;t<length_contact;t++){
      int center=centerMap[t], cofactor=rIdxMap[t];
      (*delta_buffer)(g, center)+= (*prob_weight_sum)(t,g);
      (*delta_buffer_cofactor)(g, cofactor)+=(*prob_weight_sum)(t,g);
    }
  }
  for(k=0;k<dim_features-1;k++){
    for(t=0; t<length_seq; t++){
      if (k<window_size/2 || (k>=window_size/2*2 && k<dim_features/2+window_size/2-2)){
        int idx=IdxMapCenter[t];
        if (idx<0) { continue;}
        (*feature_buffer)(k, t)= getFeatures(idx, k);
      } else {
        int idx=IdxMapCofactor[t];
        if (idx<0) continue;
        (*feature_buffer_cofactor)(k, t) = getFeatures(idx, k);
      }
    }
  }
  for ( g=0; g<num_gates; g+=SM)
  for ( k=0; k<dim_features-1; k+=SM)
  for ( t=0; t<length_seq; t+=SM)
    for ( g2=0, pdb2=&((*delta_buffer)(g,t)), 
                pdbc2=&((*delta_buffer_cofactor)(g,t)),
                pGrad2 = &m_pModel->grad[m_pModel->iFeatureStart + g*dim_features];
          g2 < (g<lastg ? SM : num_gates%SM);
          g2++, pdb2 += delta_buffer->cols,
                pdbc2 += delta_buffer_cofactor->cols,
                pGrad2 += dim_features)
    for ( k2=0, pfb2 = &((*feature_buffer)(k, t)),
                pfbc2=&((*feature_buffer_cofactor)(k, t));
          k2 < (k<lastk ? SM : (dim_features-1-dim_pairwise_features) % SM);
          k2++, pfb2 += feature_buffer->cols, 
                pfbc2 += feature_buffer_cofactor->cols)
    {
      kk=k+k2;
      for ( t2=0; t2 < (t<lasttseq ? SM : length_seq%SM); t2++) 
      {
        if (kk<window_size/2 || (kk>=window_size/2*2 && kk<dim_features/2+window_size/2-2))
          pGrad2[kk] += pdb2[t2]*pfb2[t2];
        else 
          pGrad2[kk] += pdbc2[t2]*pfbc2[t2];
      }
    }
  for ( g=0; g<num_gates; g++)
    for ( t=0, pGrad2 = &m_pModel->grad[m_pModel->iFeatureStart + g*dim_features]; 
          t<length_contact; t++)
      for ( p=0; p<dim_pairwise_features+1; p++)
        pGrad2[dim_features-dim_pairwise_features-1+p] += 
              (*prob_weight_sum)(t,g) * (*pairwise_features)(p,t);

if (proc_id==0) {
  gettimeofday(&te, 0);
  double dif = difftime (te.tv_sec,ts.tv_sec);
  //cerr << "Gradient Empirical: " << dif*1000+te.tv_usec/1000-ts.tv_usec/1000 << " miliseconds." << endl;
  gettimeofday(&ts, 0);
}
  // compute the expected values
  prob_weight_sum->Fill(0);
  for(t=0;t<length_contact;t++)
    for (l=0; l<num_states; l++)
      (*probs)(t,l) = obj_weights[t]*exp( (*arrVi)(t, l) - Partition[t] );

// out=(*gates)(t, g); prob_out = prob(t, l)*out; m_pModel->grad[l, g] -= prob_out; (*prob_weight_sum)(t,g) += prob_out*weightsL(l, g)*(1.0-out);

  for ( t=0;t<length_contact;t+=SM )
  for ( l=0; l<num_states; l+=SM )
  for ( g=0; g<num_gates; g+=SM )
    for ( t2 = 0, probs2 = &((*probs)(t,l)),
                  pPWS2 = &((*prob_weight_sum)(t,g)),
                  pGate2 = &((*gates)(t,g)); 
         t2 < (t<lastt?SM:length_contact%SM); 
         t2 ++, probs2+=probs->cols,
                pPWS2+=prob_weight_sum->cols, 
                pGate2+=gates->cols )
    {
      for ( l2 = 0, pW2=m_pModel->GetWeightLAddr(l ,g),
                    pGrad2=&m_pModel->grad[gstart + l * num_gates+g]; 
            l2 < (l<lastl?SM:num_states%SM); 
            l2++, pW2 += num_gates, pGrad2 += num_gates)
      {
        for ( g2 = 0; g2 < (g<lastg ? SM : num_gates%SM); g2++)
        {
          prob_out = probs2[l2]*pGate2[g2];
          pGrad2[g2] -= prob_out;
          pPWS2[g2] += prob_out * pW2[g2] * (1.0-pGate2[g2]);
        }
      }
    }
  delta_buffer->Fill(0);
  delta_buffer_cofactor->Fill(0);
  for(g=0;g<num_gates;g++)
  {
    for(t=0;t<length_contact;t++)
    {
      int center=centerMap[t], cofactor=rIdxMap[t];
      (*delta_buffer)(g, center) += (*prob_weight_sum)(t,g);
      (*delta_buffer_cofactor)(g, cofactor) += (*prob_weight_sum)(t,g);
    }
  }
  for ( g=0; g<num_gates; g+=SM )
  for ( k=0; k<dim_features-1; k+=SM )
  for ( t=0; t<length_seq; t+=SM )
  {
    for ( g2=0, pdb2=&((*delta_buffer)(g,t)), 
                pdbc2=&((*delta_buffer_cofactor)(g,t)), 
                pGrad2 = &m_pModel->grad[m_pModel->iFeatureStart + g*dim_features];
          g2<(g<lastg?SM:num_gates%SM); 
          g2++, pdb2+=delta_buffer->cols, 
                pdbc2+=delta_buffer_cofactor->cols, 
                pGrad2+=dim_features )
    for ( k2=0, pfb2=&((*feature_buffer)(k, t)), 
                pfbc2=&((*feature_buffer_cofactor)(k, t)); 
          k2 < (k<lastk ? SM : (dim_features - 1 - dim_pairwise_features) % SM); 
          k2++, pfb2 += feature_buffer->cols, 
                pfbc2 += feature_buffer_cofactor->cols )
    {
      kk=k+k2;
      for ( t2=0; t2<(t<lasttseq?SM:length_seq%SM); t2++) 
      {
        if (kk<window_size/2 || (kk>=window_size/2*2 && kk<dim_features/2+window_size/2-2))
          pGrad2[kk] -= pdb2[t2]*pfb2[t2];
        else 
          pGrad2[kk] -= pdbc2[t2]*pfbc2[t2];
      }
    }
  }
  for(g=0;g<num_gates;g++)
    for(t=0, pGrad2 = &m_pModel->grad[m_pModel->iFeatureStart + g*dim_features]; 
        t<length_contact; t++) {
      for (p=0; p<dim_pairwise_features+1; p++) {
        pGrad2[dim_features-dim_pairwise_features-1+p] -= 
                       (*prob_weight_sum)(t,g) * (*pairwise_features)(p,t);
      }
    }

if (proc_id==0) {
  gettimeofday(&te, 0);
  double dif = difftime (te.tv_sec,ts.tv_sec);
  //cerr << "Gradient Expected: " << dif*1000+te.tv_usec/1000-ts.tv_usec/1000 << " miliseconds." << endl;
}

  delete delta_buffer;
  delete delta_buffer_cofactor;
  delete probs,
  delete feature_buffer;
  delete feature_buffer_cofactor;
  delete prob_weight_sum;
  delete arrVi; arrVi=NULL;

  return obj;
}

void PairwiseSequence::ComputeTestAccuracy(){
  m_pModel->totalPos += length_contact;
  // comparison
  //ScoreMatrix dist_pred(length_seq, length_seq,"dist_pred");
  //dist_pred.Fill(0);
  for(int t=0; t<length_contact;t++) {
    if(obs_label[t]==predicted_label[t])
      m_pModel->totalCorrect++;
    confusion[predicted_label[t]][obs_label[t]]++;
/*
    int x=centerMap[t], y=rIdxMap[t];
    int label=predicted_label[t];
    if (label==0)
      dist_pred(x,y)=3.5;
    else if (label<m_pModel->num_labels-1)
      dist_pred(x,y)=label+3.5;
*/
  }
/*
  for (int x=m_pModel->window_size/2; x<length_seq-m_pModel->window_size/2; x++){
    for (int y=x+7; y<length_seq-m_pModel->window_size/2; y++){
      if (dist_pred(x,y)>0 && dist_pred(x,y+1)>0){
        m_pModel->totalTriangle++;
        Score diff = abs(dist_pred(x,y)-dist_pred(x,y+1));
        if (diff<4)
          m_pModel->legalTriangle++;
      }
    }
  }
*/
  //for(int t=0; t<length_contact;t++) cerr << predicted_label[t] << ":" << obs_label[t] << " "; 
  //cerr << "=="<<m_pModel->totalCorrect << endl;
}

/////////////////////  implementation of the class ContactNumberSequence ///////////////////////////////////////

DistancePotentialSequence::DistancePotentialSequence(int len, PNN1* pModel, string name) 
  : PairwiseSequence(len, pModel, name)
{
}

int DistancePotentialSequence::makeFeatures(){
  if (obs_label==NULL) 
    obs_label = new int[length_contact];
  int df = length_contact*2;
  //cerr << proc_id << "@ df=" << df << ", len_contact=" << length_contact << ", gaps_size=" << gaps.size();
  _features = new int[df];
  
  int totalContact=0;

  //cerr << " window_size=" << window_size << ", len_seq=" << length_seq << endl;
  //build features for local windows  
  int negativeLabel=0;
  //for(int t=max(s_start, window_size/2);t<min(s_end, length_seq-window_size/2);t++){
  //for(int t=window_size/2;t<length_seq-window_size/2;t++){
  for(int t=0;t<length_seq-contact_DIST;t++){
    bool bGap=false;
    for (int g=0; g<(int)gaps.size(); g++){
      if (abs(gaps[g]-t)<=window_size/2){
        bGap=true; break;
      }
    }
    if(bGap)
      continue;
    // for(int u=t+7;u<length_seq-window_size/2;u++){
    for(int u=t+contact_DIST;u<length_seq;u++){
      if (totalContact>=length_contact-1)
        break;
      bGap=false;
      for (int g=0; g<(int)gaps.size(); g++){
        if (abs(gaps[g]-u)<=window_size/2){
          bGap=true; break;
        }
      }
      if(bGap) 
        continue;
      int label = (int)(((*distances)(t,u)-3)*1);
      //int label1 = (int)(((*distances)(t,u)-7)*1), max_label=12;   // 8 states: 3-7.99, 8-8.99, ..., 13-13.99, 14-inf
      //int label1 = (int)(((*distances)(t,u)-6)*1), max_label=12; // 9 states: 3-6.99, 7-8.99, 8-8.99, ..., 13-13.99, 14-inf
      //int label1 = (int)(((*distances)(t,u)-5)*1), max_label=12; // 10 states: 3-5.99, 6-6.99, 7-7.99, 8-8.99, ..., 13-13.99, 14-inf
      //int label1 = (int)(((*distances)(t,u)-3)*1), max_label=12; // 12 states: 3-3.99, 4-4.99,5-5.99, 6-6.99, 7-7.99, 8-8.99, ..., 13-13.99, 14-inf
      int label1 = (int)(((*distances)(t,u)-4)*1), max_label=13; // 12 states: 3-4.99, 5-5.99, 6-6.99, 7-7.99, 8-8.99, ..., 13-13.99, 14-14.99, 15-inf
      //int label1 = (int)(((*distances)(t,u)-4)*1), max_label=18; // 17 states: 3-4.99, 5-5.99, 6-6.99, ..., 18-18.99, 19-19.99, 20-inf
      if (label1>0 && label<max_label)
        obs_label[totalContact]=label1;
        //obs_label[totalContact]=(label1+1)/2;
      else if (label>=0 && label<max_label)
        obs_label[totalContact]=0;
      else { 
        // when label is negative, it means distance is over 21A for 17 states 
        // or over 16A for 12 states, so it's obs_label is the last state possible
        double prob=drand48();
        if (prob>1.99999 && test_mode==0)
          continue;
        obs_label[totalContact]=num_states-1;
        negativeLabel++;
      } 
      // if (!proc_id) 
      //     cerr << "    Adding contact " << totalContact << ":" << p_name 
      //          << " (" << t << ", " << u << ") with dist=" << (*distances)(t,u)
      //          << " labeled " << obs_label[totalContact] << endl;
      centerMap[totalContact]=t;
      IdxMapCenter[t]=totalContact;
      rIdxMap[totalContact]=u;
      IdxMapCofactor[u]=totalContact;
      int pivot = totalContact*2;
      totalContact++;
      int offset1 = t-window_size/2;
      int offset2 = u-window_size/2;
      _features[pivot++]=offset1;
      _features[pivot++]=offset2;
    }
  }
  //cerr << proc_id << "@  Ratio: " << negativeLabel << "/" << totalContact << endl;
  //if (length_contact <=totalContact+2)
  //  cerr << "  ^^^ Got more contact than expected: " << length_contact << " : " << totalContact<< endl;
  delete distances; distances=NULL;
  nDataPoints=length_contact=totalContact; 
  if (length_contact==0)
     return -1; //cerr << proc_id << " Making features..." << p_name<< endl;
  predicted_label = new int[nDataPoints];
  pairwise_features = new ScoreMatrix(1, nDataPoints, "pairwise_features");
  for(int t=0;t<length_contact; t++){
    int offset1=centerMap[t], offset2=rIdxMap[t];
    (*pairwise_features)(0,t)= (offset2-offset1+0.000001)/length_seq;
  }
  obj_weights=new double[nDataPoints];
  for (int t=0; t<nDataPoints; t++)
    obj_weights[t]=1.0;

  return 0;
}

Score DistancePotentialSequence::getFeatures(int pos, int i){
  int pivot=pos*2;
  int offset1 = _features[pivot];
  int offset2 = _features[pivot+1];
  int iif=i-window_size/2*2;
  int idx = iif / m_pModel->dim_one_pos;
  Score ret=1; // when i==dim_features-3
  if (i<window_size/2*2){
    offset1+=window_size/2;
    offset2+=window_size/2;
    //if (proc_id==0) cerr << "offset1="<<offset1<<", offset2="<<offset2<<endl;
    if (offset1<window_size/2 && i==offset1)
      return 1;
    else if (length_seq-offset2==window_size/2*2-i && offset2>=length_seq-window_size/2){
      //if (proc_id==0) cerr << length_seq << "-" << offset2 << "="<<window_size/2*2<<"-" <<i << endl; 
      return 1;
    } else 
      return 0;
  } else if (idx<window_size){
    if (offset1+idx<0 || offset1+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset1+idx][iif%m_pModel->dim_one_pos];
  } else if (i<dim_features-3){
    idx-=m_pModel->window_size;
    //cerr << idx << " - " << i << "*"<<dim_features <<" : " << i%m_pModel->dim_one_pos<<" : " << length_seq << " : " << offset1 << "  : " << offset2;
    if (offset2+idx<0 || offset2+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset2+idx][iif%m_pModel->dim_one_pos];
    //cerr << " : " << ret << endl;
  } else if (i==dim_features-2)
    ret = Rg/10;
  else if (i==dim_features-1)
    ret = (offset2-offset1+0.000001)/length_seq;

  return ret;
}

/////////////////////  implementation of the class ContactNumberSequence ///////////////////////////////////////

ContactNumberSequence::ContactNumberSequence(int len, PNN1* pModel, string name) : PairwiseSequence(len, pModel, name)
{
  contact_DIST=3;
  obs_cnum=new int[len];
  memset(obs_cnum, 0, sizeof(int)*len);
}

#define _MAX_CONTACT 10
int ContactNumberSequence::makeFeatures(){
  if (obs_label==NULL) 
    obs_label = new int[length_contact];
  int df = length_contact*2;
  // cerr << proc_id << "@ df=" << df << ", len_contact=" << length_contact << ", gaps_size=" << gaps.size();
  _features = new int[df];
  
  int totalContact=0;
  int window_size = m_pModel->window_size;

  //cerr << " window_size=" << window_size << ", len_seq=" << length_seq << endl;
  int* nContacts=new int[length_seq];
  int* nContactsR=new int[length_seq];
  for(int t=0;t<length_seq-contact_DIST;t++){
      bool bGap=false;
      for (int g=0; g<(int)gaps.size(); g++){
        if (abs(gaps[g]-t)<=window_size/2){
      bGap=true; break;
    }
      }
      if(bGap)
    continue;
      memset(nContacts, 0, sizeof(int)*length_seq);
      memset(nContactsR, 0, sizeof(int)*length_seq);
    int nn=0;
      for(int u=t+contact_DIST;u<length_seq;u++){
    nContacts[u]=nContacts[u-1];
    bGap=false;
    for (int g=0; g<(int)gaps.size(); g++)
      if (abs(gaps[g]-u)<=window_size/2){ bGap=true; break; }
    if(bGap) continue;

// For predicting the contact number between residue t and all the residues from t+contact_DIST to u
    if ((*distances)(t,u)<=8.0 && (*distances)(t,u)>2) nContacts[u]++, nn++;
    if (nContacts[u]<_MAX_CONTACT-1) obs_label[totalContact]=_MAX_CONTACT-1-nContacts[u];
    else obs_label[totalContact]=0; // We merge those bins with contact number more than 9 as one bin
    //if (proc_id==0) cerr << "(" << t << ", " << u << "): " << (*distances)(t,u) << "  " << nContacts[u] << endl;
    contactStat[obs_label[totalContact]]++;

    IdxMap[t*length_seq+u]=totalContact;
    centerMap[totalContact]=t; IdxMapCenter[t]=totalContact; rIdxMap[totalContact]=u; IdxMapCofactor[u]=totalContact;
    int pivot = totalContact*2;
    totalContact++;
    int offset1 = t-m_pModel->window_size/2, offset2 = u-m_pModel->window_size/2;
    _features[pivot++]=offset1; _features[pivot++]=offset2;
      }
      for(int u=t-contact_DIST;u>=0;u--){
    nContactsR[u]=nContactsR[u+1];
    bGap=false;
    for (int g=0; g<(int)gaps.size(); g++)
      if (abs(gaps[g]-u)<=window_size/2){ bGap=true; break; }
    if(bGap) continue;
// For predicting the contact number between residue t and all the residues from t+contact_DIST to u
    if ((*distances)(t,u)<=8.0 && (*distances)(t,u)>2) nContactsR[u]++, nn++;
    if (nContactsR[u]<_MAX_CONTACT-1) obs_label[totalContact]=_MAX_CONTACT-1-nContactsR[u];
    else obs_label[totalContact]=0; // We merge those bins with contact number more than 9 as one bin
    //if (proc_id==0) cerr << "(" << t << ", " << u << "): " << (*distances)(t,u) << "  " << nContacts[u] << endl;
    contactStat[obs_label[totalContact]]++;

    IdxMap[t*length_seq+u]=totalContact;
    centerMap[totalContact]=t; IdxMapCenter[t]=totalContact; rIdxMap[totalContact]=u; IdxMapCofactor[u]=totalContact;
    int pivot = totalContact*2;
    totalContact++;
    int offset1 = t-m_pModel->window_size/2, offset2 = u-m_pModel->window_size/2;
    _features[pivot++]=offset1; _features[pivot++]=offset2;
      }
      nn=nContactsR[0]+nContacts[length_seq-1];
      obs_cnum[t]=min(_MAX_CONTACT-1, nn);
    //contactStat[nn]++;
  }
  delete nContacts; delete nContactsR; delete distances; distances=NULL;
  //if (length_contact <=totalContact+2)
  //  cerr << "  ^^^ Got more contact than expected: " << length_contact << " : " << totalContact<< endl;
  nDataPoints=length_contact=totalContact; 
  if (length_contact==0)
     return -1;
  predicted_label = new int[nDataPoints];
  pairwise_features = new ScoreMatrix(1, nDataPoints, "pairwise_features");
  for(int t=0;t<length_contact; t++){
    int offset1=centerMap[t], offset2=rIdxMap[t];
    (*pairwise_features)(0,t)= (offset2-offset1+0.000001)/length_seq;
  }
  obj_weights=new double[nDataPoints];
  for (int t=0; t<nDataPoints; t++)
    obj_weights[t]=1.0/length_seq;
  for (int t=0; t<length_seq; t++)
    obj_weights[IdxMap[t]] = obj_weights[IdxMap[length_seq+t]]
                           = obj_weights[IdxMap[(t+1)*length_seq-2]]
                           = obj_weights[IdxMap[(t+1)*length_seq-1]] 
                           = 1;
  return 0;
}

// Features: left-end vector([window_size/2]), features read from input file (dim_one_pos*window_size), right-end vector ([window_size/2]), sequence distance; sequence distance/length_seq.
Score ContactNumberSequence::getFeatures(int pos, int i){
  int pivot=pos*2;
  int offset1 = _features[pivot];
  int offset2 = _features[pivot+1];
  int iif=i-window_size/2*2;
  int idx = iif / m_pModel->dim_one_pos;
  Score ret=1; // when i==dim_features-3
  if (i<window_size/2*2){
    offset1+=window_size/2;
    offset2+=window_size/2;
    //if (proc_id==0) cerr << "offset1="<<offset1<<", offset2="<<offset2<<endl;
    if ((offset1<window_size/2 && i==offset1) || (offset2<window_size/2 && i==offset2))
      return 1;
    else if ((length_seq-offset1==window_size/2*2-i && offset1>=length_seq-window_size/2)
      || (length_seq-offset2==window_size/2*2-i && offset2>=length_seq-window_size/2)){
      //if (proc_id==0) cerr << length_seq << "-" << offset2 << "="<<window_size/2*2<<"-" <<i << endl; 
      return 1;
    } else 
      return 0;
  } else if (idx<window_size){
    if (offset1+idx<0 || offset1+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset1+idx][iif%m_pModel->dim_one_pos];
  } else if (i<dim_features-3){
    idx-=m_pModel->window_size;
    //cerr << idx << " - " << i << "*"<<dim_features <<" : " << i%m_pModel->dim_one_pos<<" : " << length_seq << " : " << offset1 << "  : " << offset2;
    if (offset2+idx<0 || offset2+idx>=length_seq) return 0; // training Secondary Structure: allow termial residues
    ret= obs_feature[offset2+idx][iif%m_pModel->dim_one_pos];
    //cerr << " : " << ret << endl;
  } else if (i==dim_features-2)
    ret = Rg/10;
  else if (i==dim_features-1)
    ret = abs(offset2-offset1+0.000001)/length_seq;

  return ret;
}

Score ContactNumberSequence::Predict()
{
  int i, j, t, startRes=1, endRes=length_seq;

  contact_DIST=3;
  ComputeGates();
  ComputeVi();
  Score obj=MAP();

  for(i=0;i<num_labels;i++) 
    for(j=0;j<num_labels;j++) 
      confusion[i][j]=0; 

  // comparison
  int* num_contact8 = new int[length_seq];
  int* pred_cnum = new int[length_seq];
  memset(num_contact8, 0, sizeof(int)*length_seq);
  memset(pred_cnum, 0, sizeof(int)*length_seq);
  int head=max(0,startRes-1), tail=min(length_seq, endRes);
  //cerr << "  --------- " << p_name << "; len_contact=" << length_contact << endl;
  //fprintf(stderr, "\n\tSeq :: (leftRes, cn)-(rightRes, cn) head=%d; tail=%d\n", head, tail);
  for (t=0; t<length_seq; t++){
    int contactH=_MAX_CONTACT-1, contactT=_MAX_CONTACT-1;
    //cerr << " ===== " << t << " :: ";
    for (i=t-contact_DIST; i>=0; i--){
      int idxH=IdxMap[t*length_seq+i];
      //cerr << idxH << ", " << predicted_label[idxH] << "," << contactH << "  ";
      if (idxH>=0 && predicted_label[idxH]<=contactH){
        contactH=predicted_label[idxH];
        //if (contactH<num_labels-1) cerr << "("<<i<<","<<num_labels-contactH-1 << ")->";
      }
    }
    //cerr << "\t||\t";
    for (i=t+contact_DIST; i<tail; i++){
      int idxT=IdxMap[t*length_seq+i];
      if (idxT>=0 && predicted_label[idxT]<=contactT){
        contactT=predicted_label[idxT];
        //if (contactT<num_labels-1) cerr << "("<<i<<","<<num_labels-contactT-1 << ")=>";
      }
    }
    //cerr << "\tpredictCN " << t << ":\t" << head << ":" << idxH << ":" << contactH << "; " << tail << ":" <<idxT << ":" << contactT << endl;
    num_contact8[t]=num_labels*2-2-(contactH+contactT);
    //cerr << "\t::" << num_contact8[t] << " vs " << obs_cnum[t] << endl;
  }
  for (t=2; t<length_seq-2; t++){
    if (num_contact8[t]==0)
      //num_contact8[t]=(num_contact8[t-3]+num_contact8[t-2]+ num_contact8[t-1]+num_contact8[t+1]+ num_contact8[t+2]+ num_contact8[t+3])/6;
      num_contact8[t]=(num_contact8[t-2]+ num_contact8[t-1]+num_contact8[t+1]+ num_contact8[t+2])/4;
      //pred_cnum[t]=(num_contact8[t-2]+ num_contact8[t-1]+num_contact8[t+1]+ num_contact8[t+2])/4;
    else if (num_contact8[t]>9)
      num_contact8[t]=9;
    pred_cnum[t]=min(_MAX_CONTACT-1,num_contact8[t]);
  }
  m_pModel->totalPos += length_seq-2;
  for(t=1; t<length_seq-1; t++) {
    if(pred_cnum[t]==obs_cnum[t])
      m_pModel->totalCorrect++;
    confusion[pred_cnum[t]][obs_cnum[t]]++;
  }

  delete num_contact8;

  return obj;
}

void _LBFGS::Report(const vector<double> &theta, int iteration, 
                    double objective, double step_length)
{
  int num_labels=m_pModel->num_labels;
  //if(iteration%5) return;
  if (!proc_id) cerr << "******* iteration=" << iteration << endl;
  m_pModel->CopyWeightsIn(theta);
  
  int tc_sum = 0, tp_sum = 0, ttT_sum = 0, tlT_sum = 0;
  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Bcast(m_pModel->weights, m_pModel->num_params, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  m_pModel->totalPos = m_pModel->totalCorrect = m_pModel->totalTriangle
                     = m_pModel->legalTriangle= 0;
  for(int i=0;i<num_labels;i++)
    for(int j=0;j<num_labels;j++) {
      confusion[i][j]=0; confusion_sum[i][j]=0;
    }
  double test_obj = 0, test_obj_sum = 0;
  for(int i=0;i<m_pModel->num_tst;i++) {
    //cerr << "Test("<< i <<")" << endl;
    test_obj+=m_pModel->testData[i]->Predict();
    delete m_pModel->testData[i]->arrVi; m_pModel->testData[i]->arrVi=NULL;
  }
//*
  test_obj = -test_obj+m_pModel->RegularizeWeights();

  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Reduce(&test_obj, &test_obj_sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Bcast(&test_obj_sum, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
// */
  double norm_w = 0;
  for(int i=0;i<theta.size();i++)
    norm_w +=theta[i]*theta[i];
  tp_sum = m_pModel->totalPos;
  tc_sum = m_pModel->totalCorrect;
  ttT_sum = m_pModel->totalTriangle;
  tlT_sum = m_pModel->legalTriangle;
  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Reduce(&m_pModel->totalPos, &tp_sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(&m_pModel->totalCorrect, &tc_sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(&m_pModel->totalTriangle, &ttT_sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(&m_pModel->legalTriangle, &tlT_sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
        MPI_Reduce(confusion, confusion_sum, 30*30, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Barrier(MPI_COMM_WORLD);
    
if(proc_id==0){  
  cout << endl << "Iteration:  " << iteration << endl;
  cout << " Weight Norm: " << sqrt(norm_w) << endl;
  cout << " Objective: " << objective << " \t test obj: \t\t" << test_obj_sum << endl;
  cout << " test ACC(MAP): " << (double) tc_sum/tp_sum <<"   " << tc_sum << "/" << tp_sum << endl;
  //cout << " test triangle: " << (double) tlT_sum/ttT_sum <<"   " << tlT_sum << "/" << ttT_sum << endl;

  cout << endl << "Confusion Matrix: " <<endl;
  memset(sum, 0, sizeof(int)*30);
  for(int i=0;i<num_labels;i++){
    double tot=0;
    for(int j=0;j<num_labels;j++){
      printf("%7d ",confusion_sum[i][j]);
      tot+=confusion_sum[i][j];
      sum[j]+=confusion_sum[i][j];
    }
    printf("%7.2f\n", 100.0*confusion_sum[i][i]/tot);
  }
  for (int i=0; i<num_labels; i++)
    printf("%7.2f ", 100.0*confusion_sum[i][i]/(0.0001+sum[i]));
  cout << endl;

  if (test_mode==1)
    return;

  string m_file = m_pModel->model_file; // model file path
  char buf[100];
  sprintf(buf,"%d",iteration);
  m_file+=buf;
  cout << m_file << endl;
  ofstream fout(m_file.c_str());
  fout << "num_params: " << m_pModel->num_params << endl;
  fout << "num_gates: " << m_pModel->num_gates << endl;
  fout << "num_labels: " << m_pModel->num_labels << endl;
  fout << "dim_one_pos: " << m_pModel->dim_one_pos << endl;
  fout << "window_size: " << m_pModel->window_size << endl;
  fout << "dim_features: " << m_pModel->dim_features << endl;
  fout << " weights: " << endl;
  for(int i=0;i<m_pModel->num_params;i++)
    fout << m_pModel->weights[i] << " ";
  fout << endl;
  fout.close();
}

  for(int i=0;i<num_labels;i++) 
    for(int j=0;j<num_labels;j++) { 
      confusion[i][j]=0; confusion_sum[i][j]=0;
    }

  m_pModel->totalPos = m_pModel->totalCorrect 
                     = m_pModel->totalTriangle = m_pModel->legalTriangle= 0;
  ttT_sum = m_pModel->totalTriangle;
  tlT_sum = m_pModel->legalTriangle;
  for(int i=0;i<m_pModel->num_data;i++){
    //cerr << proc_id << ": training data " << i << endl;
    m_pModel->trainData[i]->Predict();
  }
  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Reduce(&m_pModel->totalPos,&tp_sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(&m_pModel->totalCorrect,&tc_sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(&m_pModel->totalTriangle, &ttT_sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(&m_pModel->legalTriangle, &tlT_sum, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
        MPI_Reduce(confusion, confusion_sum, 30*30, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);

if(proc_id==0){
  cout <<" train ACC: " << (double) tc_sum/tp_sum << "   " << tc_sum << "/" << tp_sum << endl;
  //cout <<" train triangle: " << (double) tlT_sum/ttT_sum <<"   " << tlT_sum << "/" << ttT_sum << endl;
  cout << endl << "Confusion Matrix: " <<endl;
  memset(sum, 0, sizeof(int)*30);
  for(int i=0;i<num_labels;i++){
    double tot=0;
    for(int j=0;j<num_labels;j++){
      printf("%7d ",confusion_sum[i][j]);
      tot+=confusion_sum[i][j];
      sum[j]+=confusion_sum[i][j];
    }
    printf("%7.2f\n", 100.0*confusion_sum[i][i]/tot);
  }
  for (int i=0; i<num_labels; i++)
    printf("%7.2f ", 100.0*confusion_sum[i][i]/(0.0001+sum[i]));
  cout << endl;
}
}

double _LBFGS::ComputeGradient(int iteration, vector<double>&g, 
                               const vector<double> &x, bool bCalculateGate, 
                               int num_consec_small_steps)
{
  //if (test_mode==1) return;

  m_pModel->CopyWeightsIn(x);

  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Bcast(m_pModel->weights, m_pModel->num_params, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  memset(m_pModel->grad, 0, sizeof(double)*m_pModel->num_params);
  Score obj=0;
  for(int i=0;i<m_pModel->num_data;i++){
    obj+=m_pModel->trainData[i]->ComputeGradient(bCalculateGate);
  }
  obj = -obj+m_pModel->RegularizeWeights();
  //if (proc_id==0) cerr << "Gradient_Obj=" << obj << endl;
  m_pModel->RegularizeGrads();
  
  double obj_sum=0;
  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Reduce(m_pModel->grad, m_pModel->grad_sum, m_pModel->num_params, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Bcast(m_pModel->grad_sum, m_pModel->num_params, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Reduce(&obj, &obj_sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Bcast(&obj_sum, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  int nonzeroDirection=0;
  double _norm=m_pModel->CopyGradsOut(g, nonzeroDirection, num_consec_small_steps);
  if(!proc_id)
    cerr << "Norm of Gradient: " << sqrt(_norm) << "; number of steep directions: " << nonzeroDirection << endl;
  return obj_sum;
}

double _LBFGS::ComputeFunction(int iteration, const vector<double> &x)
{
  //if (test_mode==1) return 0;

  m_pModel->CopyWeightsIn(x);

  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Bcast(m_pModel->weights, m_pModel->num_params, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Barrier(MPI_COMM_WORLD);

  double obj = 0, obj_sum = 0;
  for(int i=0;i<m_pModel->num_data;i++){
    //if (proc_id==0) cerr << "Obj("<< i <<"):";
    obj += m_pModel->trainData[i]->Obj();
  }
  //cerr << "Now need to RegularizeWeights" << endl;
  obj = -obj+m_pModel->RegularizeWeights();

  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Reduce(&obj, &obj_sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Bcast(&obj_sum, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  return obj_sum;
}

void _LBFGS::Report(const string &s)
{
  if(!proc_id) cerr << s << endl;
}

void PNN1::SetSeed()
{
  unsigned int randomSeed=0;
  ifstream in("/dev/urandom",ios::in);
  in.read((char*)&randomSeed, sizeof(unsigned)/sizeof(char));
  in.close();
          
  unsigned id=getpid();
  randomSeed=randomSeed*randomSeed+id*id;

  //we can set the random seed at only the main function
  srand48(randomSeed);
  srand(randomSeed);
}

int Format;
double regularizer;
int init_from_file;
string init_file;
int num_test_set;
string prefix_model_file;
string model_path;
void PNN1::SetParameters(int model, int w_size, int n_labels, int n_gates, 
                         int n_local, int n_pairwise)
{
  this->model = (MODEL) model;
  SetSeed();
  window_size = w_size; //best 13
  num_labels = n_labels;
  num_states = num_labels; // S1,...,Sk
  num_gates = n_gates; //best 20
  dim_one_pos = n_local;  // dim of local feature + 1 (for bias)
  dim_pairwise_features = n_pairwise;
  if (this->model==singleton) // Secondary Structure
    dim_features = window_size/2*2 + window_size*dim_one_pos + 1; 
  else if (this->model==pairwise) // constance 1 
    dim_features = window_size/2*2 + (2*window_size*dim_one_pos) + n_pairwise + 1; 
  else // 2+1 means Rg, chain distance(cd)/length_seq, and constance 1 
    dim_features = window_size/2*2 + (2*window_size*dim_one_pos) + 2 + 1; 
  
  if (proc_id==0)
    cerr << " dim_features =" << dim_features << "=" << window_size/2 << "*" << 2 
         << "+" << window_size << "*" <<dim_one_pos<< "+" << 1 << endl;
  iGateStart = 0;
  iFeatureStart = iGateStart+num_states*num_gates;
  num_params = iFeatureStart + dim_features*num_gates;

  weights = new double[num_params];
  memset(weights, 0, sizeof(double)*num_params);
  grad = new double[num_params];  
  memset(grad, 0, sizeof(double)*num_params);
  grad_sum = new double[num_params];
  reg = new double[num_params];  
  memset(reg, 0, sizeof(double)*num_params);
  double r=regularizer/num_procs; // regularization coefficients
  int i;
  for(i=0;i<iGateStart;i++) 
    reg[i] = r*1;
  for(;i<iFeatureStart;i++) 
    reg[i] = r*1;
  for(;i<num_params;i++)
    reg[i] = r;
}

void PNN1::convertModel(string oldModel, string newModel)
{
  ifstream fin(oldModel.c_str()); // no file header
  char buf[10240];
  string tmp;
  int data;
  for (int i=0; i<8; i++) {
    if (!fin.getline(buf,10240))
      break;
    else {
      istringstream si(buf);
      si >> tmp >> data;
      if (tmp=="num_params:")
        num_params = data;
      else if (tmp=="num_gates:")
        num_gates = data;
      else if (tmp=="window_size:")
        window_size = data;
      else if (tmp=="dim_features:")
        dim_features = data;
      else if (tmp=="num_labels:")
        num_labels = data;
      else if (tmp=="dim_one_pos:")
        dim_one_pos = data;
      cerr << tmp << "::::"<<data<<endl;
    }
  }
  weights=new double[num_params];
  for(int i=0;i<num_params;i++)
    fin >> weights[i];
  fin.close();

  int dfeatures = window_size/2*2+(2*window_size*dim_one_pos)+2+1;
  iFeatureStart = num_labels*num_gates;
  int nparams = iFeatureStart + dfeatures*num_gates;

  ofstream fout(newModel.c_str());
  fout << "num_params: " << nparams << endl;
  fout << "num_gates: " << num_gates << endl;
  fout << "num_labels: " << num_labels << endl;
  fout << "dim_one_pos: " << dim_one_pos << endl;
  fout << "window_size: " << window_size << endl;
  fout << "dim_features: " << dfeatures << endl;
  fout << " weights: " << endl;

  for(int i=0;i<iFeatureStart;i++)
    fout << weights[i] << " ";
  for(int i=0; i<num_gates; i++){
    for (int j=0; j<window_size/2*2; j++)
      fout << 0.0 << " ";
    int ifs = iFeatureStart+i*dim_features;
    for (int j=0; j<dim_features; j++){
      int idx=ifs+j;
      fout << weights[idx] << " ";
    }
    fout << 0.0 << " ";
  }
  fout << endl;

  fout.close();
}

// return the total number of valid columns
int PNN1::getColumns(string columns, int* column_flags)
{
  memset(column_flags, 0, sizeof(int)*MAX_COL);
  if (columns.length()>0) {
    //cerr << "|" << columns << "|" << endl;
    istringstream si(columns);
    char c=' ';
    int col1,col2;
    while (!si.eof()){
      si >> col1 >> c;
      //cerr << "**" << col1 << ":" << c << endl;
      if (c==','){
        if (col1>=0 && col1<MAX_COL) column_flags[col1]=1;
      } else if (c=='-'){
        si >> col2 >> c;
        if (col1>=0 && col1<MAX_COL){
          for (int i=col1; i<=min(MAX_COL-1, col2); i++)
            column_flags[i]=1;
        }
      } else {
        if (col1>=0 && col1<MAX_COL) column_flags[col1]=1;
        break;
      }
    }
    //cerr << "selected columns: ";
    //for (int i=0; i<MAX_COL; i++)
    //  if (column_flag[i]==1) cerr << i << ",";
    //cerr << endl;
  }
  int total_col=0;
  for (int i=0; i<MAX_COL; i++)
    total_col+=column_flags[i];
  return total_col;
}

void PNN1::Initialize(int model, string model_dir, int w_size, int n_labels, int n_gates, string columns, string pairwise_columns, string input_f, string split_file)
{
  char fname[100];
  if(init_from_file)
    sprintf(fname, "%s/%s_model.r%f_g%i_w%i_f_", model_path.c_str(), prefix_model_file.c_str(), regularizer, n_gates, w_size);
  else
    sprintf(fname, "%s/%s_model.r%f_g%i_w%i_", model_path.c_str(), prefix_model_file.c_str(), regularizer, n_gates, w_size);
  model_file = model_dir+fname;
  int n_local = 1 + getColumns(columns, column_flags);
  int n_pairwise = getColumns(pairwise_columns, pairwise_column_flags);
  SetParameters(model, w_size, n_labels, n_gates, n_local, n_pairwise);
  trainData.clear();
  testData.clear();
  //cerr << "Format is " << Format << endl;
  if (Format==0)
    LoadData(input_f, split_file);
  else
    LoadDataSVMFormat(input_f);

if(!proc_id){  
  cerr <<"num_params = " << num_params << endl;
  if(init_from_file){
    cerr << "  init_file: " << init_file << endl;
    ifstream fin(init_file.c_str()); // no file header
    char buf[10240];
    string tmp;
    int data;
    for (int i=0; i<7; i++) {
      fin.getline(buf,10240);
      if (fin.gcount()==0)
        break;
      else {
        istringstream si(buf);
        si >> tmp >> data;
        if (tmp=="num_params:")
          num_params = data;
        else if (tmp=="num_gates:")
          num_gates = data;
        else if (tmp=="window_size:")
          window_size = data;
        else if (tmp=="dim_features:")
          dim_features = data;
        else if (tmp=="num_labels:")
          num_labels = data;
        else if (tmp=="dim_one_pos:")
          dim_one_pos = data;
        cerr << tmp << "::::"<<data<<endl;
      }
    }
    cerr << strlen(buf) << "  Last 10 weights are : " << endl << "    ";
    for(int i=0;i<num_params;i++) {
      fin >> weights[i];
      if (i>num_params-10) cerr << weights[i] << ", ";
    }
    fin.close();
    cerr << endl;
/*
    cerr << "###&&&&*** Weights for features start at #" << iFeatureStart << endl;
    for (int g=0; g<num_gates; g+=5){
      for (int i=0; i<window_size+3; i++)
        cerr << GetWeightG(g,i) << ", ";
      cerr << " .... " << GetWeightG(g,dim_features-3) << ", " 
           << GetWeightG(g,dim_features-1) << ", " << GetWeightG(g,dim_features-1) << endl;
    }
    exit(0);
// */
  } else {
    cerr <<"nnz_param = " << nnz_param << endl;
    vector<double> w(nnz_param,0);
    for (int i=0; i<nnz_param; i++){
      w[i]=(drand48()-drand48())*0.1;
      //cerr << i << ": " << w[i] << endl;
    }
    CopyWeightsIn(w);
  }
}
  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Bcast(weights, num_params, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Bcast(reg, num_params, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  if (!proc_id) cerr << "Making features..." << endl;
  memset(contactStat, 0, sizeof(int)*30);
  memset(contactStat_sum, 0, sizeof(int)*30);
  for(int i=0;i<num_data;i++){
    //if (!proc_id) cerr << "Making training features..." << i << trainData[i]->p_name<< endl;
    if (trainData[i]->makeFeatures()==-1){
      trainData.erase(trainData.begin()+i);
      i--, num_data--;
    }
  }
  MPI_Barrier(MPI_COMM_WORLD);
        MPI_Reduce(contactStat, contactStat_sum, 30, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Barrier(MPI_COMM_WORLD);
  if (proc_id==0 && model==contact_number){
    //for (int i=0; i<num_labels; i++)
    //  cerr << i << ": " << contactStat_sum[i] << endl;
  }
  memset(contactStat, 0, sizeof(int)*30);
  memset(contactStat_sum, 0, sizeof(int)*30);
  for(int i=0;i<num_tst;i++){
    //if (!proc_id) cerr << "Making test features..." << i << testData[i]->p_name<< endl;
    if (testData[i]->makeFeatures()==-1){
      testData.erase(testData.begin()+i);
      i--, num_tst--;
    }
  }
  MPI_Barrier(MPI_COMM_WORLD);
        MPI_Reduce(contactStat, contactStat_sum, 30, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Barrier(MPI_COMM_WORLD);
/*
  if (proc_id==0 && model==contact_number){
    cerr << "\ttest Data" << endl;
    for (int i=0; i<num_labels; i++)
      cerr << i << ": " << contactStat_sum[i] << endl;
    //exit(0);
  }
  if (proc_id==1){
    cerr << "  First 10 weights are : " << endl << "    ";
    for (int i=0; i<10; i++)
      cerr << weights[i] << ", ";
    cerr << endl;
    cerr << "  Last 10 weights are : " << endl << "    ";
    for (int i=num_params-11; i<num_params; i++)
      cerr << weights[i] << ", ";
    cerr << endl;
  }
*/
}

void myunexpected () {
  cerr << "unexpected handler called " << proc_id << endl;
  throw;
}

void PNN1::LoadData(string input, string split_file)
{
  vector<int> myseq_train, myseq_test;
  if (split_file.length()){
    ifstream sf;
    string line;
    sf.open(split_file.c_str(), ifstream::in);
    for(int i=0;i<=proc_id;i++)
      getline(sf, line);
    int lastspace=0, i=0, bTest=false;
    //if (proc_id==0) cerr << "=" << line << endl;
    for (; i<line.size(); i++){
      if (line[i]==';'){
        bTest=true;
        lastspace = i+1;
        continue;
      }
      if (line[i]==' '){
        string token=line.substr(lastspace, i-lastspace);
        //if (proc_id==1) cerr << "token=" << token << endl;
        if (bTest)
          myseq_test.push_back(atoi(token.c_str()));
        else
          myseq_train.push_back(atoi(token.c_str()));
        lastspace = i+1;
      }
    }
    string token=line.substr(lastspace, i-lastspace);
    //if (proc_id==1) cerr << "token=" << token << endl;
    if (token.size())
      myseq_test.push_back(atoi(token.c_str()));
  } 

  // Read the training data
  ifstream trn_in(input.c_str());
  trn_in >> num_data;
  if(!proc_id){
    cerr << input << ": gateStart="<<iGateStart<<"; ifs="<<iFeatureStart 
         << " - Load Data: "<< dim_one_pos << ";" << num_data 
         << "; model=" << model << "; pairwise=" << dim_pairwise_features << endl;
  }
  //num_data = 200;
  int total_data = num_data;
  vector<SEQUENCE*> DATA;

  bool title_line_read=false;
  double dist, tmp;
  int a1, a2, length_s, length_contact, real_dim_one_pos, gap;
  string stmp, sline;
  SEQUENCE *seq;
  getline(trn_in, sline);
  for(int i=0;i<min(16000,num_data);i++){
    bool exist_train=false, exist_test=false;
    seq=NULL;
    if (split_file.length()){
      for (int k=0; k<myseq_train.size(); k++){
        if (i==myseq_train[k]){
          exist_train=true;
          break;
        }
      }
      for (int k=0; k<myseq_test.size(); k++){
        if (i==myseq_test[k]){
          exist_test=true;
          break;
        }
      }
    } else if(i%num_procs == proc_id || test_mode==1) 
      exist_train = true;

    if (!title_line_read)
      getline(trn_in, sline);
    else
      title_line_read = false;
    // construct a new sequence
    //if(!proc_id) cerr << "sequence title: " << sline << endl;
    switch (model)
    {
      case singleton:
      { // e.g. Secondary Structure
        istringstream iss(sline);
        iss >> length_s >> stmp; 
        //getline(trn_in, stmp);
        if (exist_train || exist_test)
          seq = new SingletonSequence(length_s, this, stmp);
        break;
      }
      case pairwise:
      {
        istringstream iss(sline);
        iss >> length_s >> length_contact >> stmp;
        //getline(trn_in, stmp);
        if (exist_train || exist_test)
          seq = new PairwiseSequence(length_s, this, stmp);
        break;
      } 
      case DistancePotential:
      {
        istringstream iss(sline);
        iss >> length_s >> length_contact >> stmp;
        //getline(trn_in, stmp);
        if (exist_train || exist_test)
          seq = new DistancePotentialSequence(length_s, this, stmp);
        break;
      } 
      case contact_number:
      {
        istringstream iss(sline);
        iss >> length_s >> length_contact >> stmp;
        //getline(trn_in, stmp);
        if (exist_train || exist_test)
          seq = new ContactNumberSequence(length_s, this, stmp);
        break;
      } 
    }
    // if(!proc_id) cerr << proc_id << " *** " << i<<":"<<length_s << ", Contact:" << length_contact << " " << stmp <<endl;

    for(int j=0;j<length_s;j++){
      getline(trn_in, sline); //cerr << "Contact " << j << ":";
      if (!exist_train && !exist_test)
        continue;
      istringstream iss(sline);
      int d=0;
      for (int k=0; !iss.eof(); k++) {
        iss >> tmp;
        if (column_flags[k]==1) {
          seq->obs_feature[j][d++] = tmp; //cerr << tmp <<" ";
        }
      } //cerr << endl;
      real_dim_one_pos=d+1;
      seq->obs_feature[j][d] = 1;
    }
    //if(!proc_id) cerr << proc_id << ": Data of residue " << i <<endl;
    
    if (model==singleton){
      for(int j=0;j<length_s;j++){
        getline(trn_in, stmp); //cerr << "Contact " << j << ":";
        if (!exist_train && !exist_test)
          continue;
        if (stmp[0]=='#'){
          gap=atoi(stmp.substr(5).c_str());
          seq->gaps.push_back(gap);
        } else {
          istringstream iss(stmp);
          iss >> a1; // if (!proc_id) cerr << a1 << "(" << j << ")" << endl;
          seq->obs_label[j]=a1;
        }
      } //exit(0);
    } else { // pairwise
      getline(trn_in, stmp);
      for(int j=0;j<length_contact;j++){
        getline(trn_in, stmp); //cerr << "Contact " << j << ":";
        if (*(stmp.end()-1)=='#') {
          sline = stmp;
          title_line_read=true;
          break;
        }
        if (!exist_train && !exist_test)
          continue;
        if (stmp[0]=='#'){
          gap=atoi(stmp.substr(5).c_str());
          seq->gaps.push_back(gap);
        } else if (stmp.size()>0){
          istringstream iss(stmp);
          iss >> dist >> a1 >> a2; 
          try {
            seq->distances->layers=proc_id;
            seq->distances->id=i;
            seq->distances->helper=dist;
            if (a1<length_s && a2<length_s){
              (*seq->distances)(a1,a2)=dist;
              (*seq->distances)(a2,a1)=dist;
            } else 
              cerr << proc_id << " ::" << stmp << ":: " << seq->p_name 
                   << " dist" << "(" << a1 << "," << a2 << ")=" << dist << endl;
            //if (!proc_id && exist_test)
            //  cerr << (*seq->distances)(a1,a2) << "(" << a1 << "," << a2 << ") " << seq->p_name << endl;
            for (int k=0, d=0; !iss.eof(); k++) {
              iss >> tmp;
              if (pairwise_column_flags[k]==1) {
                (*seq->obs_pairwise_features)(a1, a2, d++) = tmp; //cerr << tmp <<" ";
              }
            } //cerr << endl;
          } catch (bad_exception e){
            cerr << proc_id << " " << seq->p_name << " dist" << "(" << a1 << "," << a2 << ")=" << dist << endl;
            exit(0);
          }
        }
      }
    }
    if (split_file.length()){
      if (exist_train && length_s>=40 && length_s<=650){
        trainData.push_back(seq);
        if (proc_id==1) cerr << " Inserting training " << i << " length=" << seq->length_contact << endl;
      } else if (exist_test && length_s>=40 && length_s<=650){
        if (proc_id==1) cerr << " Inserting testing " << i << " length=" << seq->length_contact << endl;
        testData.push_back(seq);
      } else{
        delete seq; continue;
      }
    } else {
      //if(!proc_id) cerr << proc_id << " Inserting sequence... " << i <<endl;
      if(i%num_procs != proc_id && test_mode==0) { 
        //cerr << proc_id << "deleting..." << seq << ". ";
        delete seq; 
        continue;
      } else if (length_s>=40 && length_s<=650){
        if (proc_id==6) cerr << " Inserting " << i << " length=" << seq->length_seq << endl;
        DATA.push_back(seq); 
      }
    }
    //cerr << proc_id << " done with " << i << endl;
  }
  dim_one_pos = real_dim_one_pos;
  trn_in.close();
  nnz_param=num_states*num_gates + dim_features*num_gates;
  if(!proc_id) {
    cerr << "num_gates=" << num_gates << endl;
    cerr << "num_params=" << num_params << "; nnz_param=" << nnz_param
         << "; num_states=" << num_states << "; dim_features=" << dim_features << endl;
  }

  num_data = DATA.size();
if(!proc_id){
    cerr << "total_data="<<total_data<< endl;
  }
  int ltest=0, ltrain=0;
  if (split_file.length()){
    for (int i=0; i<trainData.size(); i++)
      ltrain+=trainData[i]->length_contact;
    for (int i=0; i<testData.size(); i++)
      ltest+=testData[i]->length_contact;
  } else {
    for(int i=0, idx=0;i<num_data;i+=1,idx++){
      if (test_mode==1)
        testData.push_back(DATA[i]);
      else if (num_test_set>0){
        if ( i >= (total_data-num_test_set) / num_procs){
          testData.push_back(DATA[i]); 
          ltest += testData[testData.size()-1]->length_contact;
        } else {
          trainData.push_back(DATA[i]);
          ltrain+=trainData[trainData.size()-1]->length_contact;
        }
      } else if (idx%4==3){
        testData.push_back(DATA[i]);
        ltest += testData[testData.size()-1]->length_contact;
      } else {
        trainData.push_back(DATA[i]);
        ltrain += trainData[trainData.size()-1]->length_contact;
      }
      //cerr << i << " " << DATA[i]->length_seq << endl;
    }
  }
  num_tst=testData.size();
  num_data=trainData.size();
  if (proc_id==206){
    cerr << "DataForProc " << proc_id << ": #data="<< num_data 
         << "(" << ltrain << "); #tst=" << num_tst << "(" << ltest << ")" << endl;
  }
  MPI_Barrier(MPI_COMM_WORLD);
}

void PNN1::LoadDataSVMFormat(string input)
{
  if (model!=singleton)
    return; // svm format only works for singletons

  if (proc_id==0) cerr << " Read from " << input << endl;
  // Read the training data
  int total_data = 0;
  string sline, token;
  vector<string> buffer;
  ifstream trn_in(input.c_str());
  while (!trn_in.eof()) {
    getline(trn_in, sline);
    buffer.push_back(sline);
    total_data ++;
  }
  num_data = num_procs*64;
  vector<SEQUENCE*> DATA;
  double tmp;
  int length_s=total_data/num_data, real_dim_one_pos;
  if(!proc_id){
    cerr << proc_id<<":"<<input<<": gateStart="<<iGateStart<<"; ifs="<<iFeatureStart
         << " - Load Data: "<<dim_one_pos<<";"<<num_data<<"; model="<<model<<"; "
         << total_data << endl;
  }
  SEQUENCE *seq;
  for(int i=0;i<num_data;i++){
    if(i%num_procs != proc_id && test_mode==0) 
      continue;

    // construct a new sequence
    seq = new SingletonSequence(length_s, this, "");
    for(int j=0;j<length_s;j++){
      sline = buffer[i*length_s+j];
      int lastspace=0, lastcolon=0, idx=0, k=0, labeled=0;
      for ( ; k<sline.size(); k++) {
        if (sline[k]==' '){
          if (!labeled) {
            string token=sline.substr(lastspace, k-lastspace);
            seq->obs_label[j]=atoi(token.c_str())-1;
            //if (j==0 && proc_id==0) cerr << endl << seq->obs_label[j] <<" ";
            labeled=1;
          } else {
            if (lastcolon==0){
              idx++;
              token=sline.substr(lastspace, k-lastspace);
            } else 
              token=sline.substr(lastcolon, k-lastcolon);
            seq->obs_feature[j][idx-1] = atof(token.c_str()); 
            //if (j==0 && proc_id==0) cerr << idx << ":" << seq->obs_feature[j][idx-1] <<" ";
          }
          lastspace = k+1;
        } else if (sline[k]==':') {
          string token=sline.substr(lastspace, k-lastspace);
          idx = atoi(token.c_str());
          lastcolon = k+1;
        }
      }
      if (lastcolon==0){
        idx++;
        token=sline.substr(lastspace, k-lastspace);
      } else 
        token=sline.substr(lastcolon, k-lastcolon);
      seq->obs_feature[j][idx-1] = atof(token.c_str()); 
      //if (j==0 && proc_id==0) cerr << idx << ":" << seq->obs_feature[j][idx-1] <<" ";
      real_dim_one_pos=idx+1;
      seq->obs_feature[j][idx] = 1;
    }
    if (proc_id==0) 
      cerr << endl << " Inserting " << i << " length=" << seq->length_seq << endl;
    DATA.push_back(seq); 
    //cerr << proc_id << " done with " << i << endl;
  }
  dim_one_pos = real_dim_one_pos;
  trn_in.close();
  nnz_param=num_states*num_gates + dim_features*num_gates;
  if(!proc_id) {
    cerr << "num_gates=" << num_gates << endl;
    cerr <<"num_params=" << num_params << "; nnz_param="<<nnz_param
         <<"; num_states=" << num_states << "; dim_features=" << dim_features<<endl;
  }

  num_data = DATA.size();
if(!proc_id){
    cerr << "total_data="<<total_data<< endl;
  }
  int ltest=0, ltrain=0;
  for(int i=0, idx=0;i<num_data;i+=1,idx++){
    if (test_mode==1)
      testData.push_back(DATA[i]);
    else if (num_test_set>0){
      if (i>=(total_data-num_test_set)/num_procs){
        testData.push_back(DATA[i]); ltest+=testData[testData.size()-1]->length_seq;
      } else {
        trainData.push_back(DATA[i]); ltrain+=trainData[trainData.size()-1]->length_seq;
      }
    } else if (idx%4==3){
      testData.push_back(DATA[i]); ltest+=testData[testData.size()-1]->length_seq;
    } else {
      trainData.push_back(DATA[i]); ltrain+=trainData[trainData.size()-1]->length_seq;
    }
    //cerr << i << " " << DATA[i]->length_seq << endl;
  }
  
  num_tst=testData.size();
  num_data=trainData.size();
  if (proc_id==206) {
    cerr << "DataForProc " << proc_id << ": #data="<< num_data 
         << "(" << ltrain << "); #tst=" << num_tst << "(" << ltest << ")" << endl;
  }
  MPI_Barrier(MPI_COMM_WORLD);
}

void PNN1::CopyWeightsOut(vector<double>& w)
{
  for (int l=0; l<nnz_param; l++)
    w[l] = weights[l];
}

void PNN1::CopyWeightsIn(const vector<double>& w)
{
  for (int l=0; l<nnz_param; l++)
    weights[l] = w[l];
}

double PNN1::CopyGradsOut(vector<double>& g, int& nonzeroDirection, int num_consec_small_steps)
{
  nonzeroDirection=0;
  double _norm=0;
  for(int l=0; l<nnz_param;l++){ 
    g[l]=grad_sum[l]; 
    //if (g[l]!=0)
    if (abs(g[l])>100){
      if (num_consec_small_steps>0){
        g[l]=0;
      }
      // if (proc_id==0 && l>30000 && l<100000) cerr << "Grad: " << g[l] << " " << l << endl;
      nonzeroDirection++;
    }
    _norm +=g[l]*g[l]; 
  }

  return _norm;
}

double PNN1::RegularizeWeights()
{
  double r=0;
  for (int l=0; l<nnz_param; l++)
    r+=weights[l]*weights[l]*reg[l];
  return r;
}

void PNN1::RegularizeGrads()
{
  for (int l=0;l<nnz_param; l++)
    grad[l]=-grad[l]+weights[l]*reg[l]*2;
}

void Usage()
{
  cerr << "PNN1 [-i input_file] [-d data_dir] [-w windowSize] [-s num_labels] "
       << "[-g gates] [-m model] [-l columncs] [-p pairwise_features] "
       << "[-r regularizer] [-F format] [-f init_model_file] [-t num_test] "
       << "[-S split_file] [-M model_path] [-R newModelFile] [-P model_file_prefix]\n";
  cerr << "\t-m model: singleton, pairwise, DistancePotential or contact_number"
       << endl;
  cerr << "\t-l columns: select those columns to calculate in the singleton "
       << "features section. Note: the columns start from 0. For example, "
       << "-l0,3-6 means columns 0, 3, 4, 5 and 6." << endl;
  cerr << "\t-p pairwise_feaures : select those columns to calculate in the "
       << "pairwise feature section. Similar to -l option. Default value is -1, meaning no feature." << endl;
  cerr << "\t-F format: 0 denotes sequence format; 1 denotes svm format. Default is 0." << endl;
  cerr << "\t-t num_test: number of sequences used as test set. It's counted "
       << "from the end of the input file. If -t option is not used, the 3rd "
       << "of every 4 sequences is included in the test set by default." << endl;
  cerr << "  Example: " << endl;
  cerr << "     mpirun -np 4 ./PNN1 -i ss70 -d . -w 5 -s 3 -g 20 -l 0-39 "
       << "-r 0.001 -m singleton -t 100" << endl;
  cerr << "  This command uses 4 cores to run a 3 state secondary structure "
       << "classfier. " << endl;
  cerr << "  The window size is set to be 5, using 20 gates, 40 features for "
       << "each residue, i.e. column 0 to 39. The regularizer is set to be 0.001." << endl;
  cerr << "  ss70 file has 513 proteins sequences. We use 100 of them from the "
       << "end as the test set." << endl;
}

int main(int argc, char **argv){

  MPI_Init(&argc, &argv);

  // the command line must be:
  if (argc < 3) {
    Usage();
    return 0;
  }

  Format=0;
  num_test_set=0;
  regularizer=0.5;
  init_from_file=0;
  prefix_model_file="";
  model_path="model";
  extern char* optarg;
  char c=0;
  int w_size = 9; //best 13
  int n_labels = 101, n_gates = 20; //best 20
  int model=1;
  string model_dir="./", split_file="", input_f="ss70", newModelFile="", 
         columns="", pairwise_columns="";
  while((c=getopt(argc,argv,":w:s:g:d:i:l:p:r:F:f:t:m:M:S:P:T:R:"))!=EOF){
    switch(c)
    {
    case 'i':
      input_f=optarg;
      break;
    case 'w':
      w_size=atoi(optarg);
      break;
    case 's':
      n_labels=atoi(optarg);
      break;
    case 'g':
      n_gates = atoi(optarg);
      break;
    case 'l':
      columns=optarg;
      break;
    case 'p':
      pairwise_columns=optarg;
      break;
    case 'd':
      model_dir = optarg;
      if (model_dir[model_dir.size()-1] != '/')
        model_dir += "/";
      //cerr << "model_dir=" << model_dir << endl;
      break;
    case 'F':
      Format = atoi(optarg);
      //cerr << " Format=" << Format << "=" << optarg << endl;
      break;
    case 'r':
      regularizer=atof(optarg);
      break;
    case 'f':
      init_from_file=1;
      init_file = optarg;
      break;
    case 't':
      num_test_set=atoi(optarg);
      break;
    case 'm':
      if (strcasecmp(optarg, "singleton")==0)
        model = 0;
      else if (strcasecmp(optarg, "pairwise")==0)
        model = 1;
      else if (strcasecmp(optarg, "DistancePotential")==0)
        model = 2;
      else if (strcasecmp(optarg, "contact_number")==0)
        model = 3;
      break;
    case 'M':
      model_path = optarg;
      break;
    case 'S':
      split_file = optarg;
      break;
    case 'P':
      prefix_model_file=optarg;
      break;
    case 'T':
      test_mode=atoi(optarg);
      break;
    case 'R':
      newModelFile=optarg;
      break;
    default:
      Usage();
      exit(-1);
    }
  }

  //init_from_file=1;
  //init_file = "modelRg3w15";
  //init_from_file=0;
  //init_file = "";
  //input_f="trainDistance5247HMM"; // special case for pseudoIsing/PIcnRg1_w15.pbs
  //input_f="trainDistance5247CIMH"; // special case for pseudoIsing/PIcnRg1_w15.pbs
  //n_labels=12; // special case for pseudoIsing/PIcnRg1_w15.pbs
  //n_labels=17; // special case for pseudoIsing/PIcnRg1_w15.pbs

  PNN1 cnfModel;
  if (newModelFile!=""){
    cnfModel.convertModel(init_file, newModelFile);
    exit(0);
  }
  MPI_Comm_rank(MPI_COMM_WORLD, &proc_id);
  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);

  set_unexpected (myunexpected);

  cnfModel.Initialize(model, model_dir, w_size, n_labels, n_gates, columns,
                      pairwise_columns, input_f, split_file);
  if(!proc_id){
    cerr << "Initialization Finished!" << endl;
    cerr << "num_data = " << cnfModel.num_data*num_procs << endl;
    cerr << "num_tst = " << cnfModel.num_tst*num_procs << endl;
  }

  _LBFGS* lbfgs = new _LBFGS(&cnfModel);
  vector<double> w_init(cnfModel.nnz_param,0);
  cnfModel.CopyWeightsOut(w_init);
  //for (int i=0; i<20; i++) cerr << w_init[i]<<","; cerr << endl;//exit(0);
  MPI_Barrier(MPI_COMM_WORLD);
  int max_iter=8000;
  if (test_mode==1){
    lbfgs->Report(w_init, 0, 0, 0);
  } else{
    lbfgs->LBFGS(w_init, max_iter);
    //lbfgs->TrainOnline(w_init, max_iter);
  }
  MPI_Barrier(MPI_COMM_WORLD);
  for(int i=0;i<cnfModel.num_data;i++)
    delete cnfModel.trainData[i];
  for(int i=0;i<cnfModel.num_tst;i++)
    delete cnfModel.testData[i];
  MPI_Finalize();
  return 0;
}
